## Callin Switzer
## 17 Nov 2017
## Multilevel model to visualize bees'
## behavior on the artificial pollen system

#install packages
ipak <- function(pkg){
     new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
     if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
     sapply(pkg, require, character.only = TRUE)
}

packages <- c("ggplot2", "reshape2", 'lme4', 'sjPlot', "multcomp", "plyr", "effects", "viridis")
ipak(packages)
## Loading required package: ggplot2
## Loading required package: reshape2
## Loading required package: lme4
## Loading required package: Matrix
## Loading required package: sjPlot
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.2.10
## Current Matrix version is 1.2.11
## Please re-install 'TMB' from source or restore original 'Matrix' package
## Loading required package: multcomp
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
## Loading required package: plyr
## Loading required package: effects
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
## Loading required package: viridis
## Loading required package: viridisLite
##  ggplot2 reshape2     lme4   sjPlot multcomp     plyr  effects  viridis 
##     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE
# set ggplot theme
theme_set(theme_classic() + theme(axis.text=element_text(colour="black")))

# define data and figure directories
dataDir <- "/Users/cswitzer/Dropbox/SonicationBehavior/SonBehData"
figDir <- "/Users/cswitzer/Dropbox/SonicationBehavior/SonBehFigs"

print(paste("last run ", Sys.time()))
## Warning in as.POSIXlt.POSIXct(x, tz): unknown timezone 'zone/tz/2017c.1.0/
## zoneinfo/America/Los_Angeles'
## [1] "last run  2018-01-11 17:39:46"
print(R.version)
##                _                           
## platform       x86_64-apple-darwin15.6.0   
## arch           x86_64                      
## os             darwin15.6.0                
## system         x86_64, darwin15.6.0        
## status                                     
## major          3                           
## minor          4.2                         
## year           2017                        
## month          09                          
## day            28                          
## svn rev        73368                       
## language       R                           
## version.string R version 3.4.2 (2017-09-28)
## nickname       Short Summer

Read in data and double check it

sl <- read.csv(file.path(dataDir, '01_CombinedTrials_cleaned.csv'))

# make sure hive is a factor
sl$hive = as.factor(sl$hive)
unique(sl$hive)
## [1] 4 3 5
## Levels: 3 4 5
# number of trials per bee
xtabs(~sl$beeCol)
## sl$beeCol
##             blue             gold          goldred            green 
##               36               54                3               13 
##             lime         limeblue         limegold        limegreen 
##               28              530               54               50 
##       limeorange         limepink       limepurple limepurpleyellow 
##               84               51              352               58 
##          limered       limesilver        limewhite       limeyellow 
##             3474                3              101               52 
##           orange       orangeblue      orangegreen       orangepink 
##              501               50               50               50 
##     orangepurple           purple          redblue         redgreen 
##               50                9              673              530 
##          redpink        redpurple           silver            white 
##              677              218               26              283 
##        whiteblue        whitegold       whitegreen      whiteorange 
##             2786             1225               39             1122 
##        whitepink      whitepurple         whitered      whiteyellow 
##             1184               55             1344              157 
##       yellowblue      yellowgreen     yelloworange       yellowpink 
##             1587              532             2113             2309 
##     yellowpurple        yellowred 
##             1243              547
hist(xtabs(~sl$beeCol))

max(xtabs(~sl$beeCol))
## [1] 3474
mean(xtabs(~sl$beeCol))
## [1] 578.6429
sd(xtabs(~sl$beeCol))
## [1] 829.7747
# number of visits for first trial
xtabs(~sl$beeCol[sl$trialNum == 1])
## sl$beeCol[sl$trialNum == 1]
##             blue             gold          goldred            green 
##               36               54                3                9 
##             lime         limeblue         limegold        limegreen 
##               28               50               54               50 
##       limeorange         limepink       limepurple limepurpleyellow 
##               33               51               53               58 
##          limered       limesilver        limewhite       limeyellow 
##               50                3               50               52 
##           orange       orangeblue      orangegreen       orangepink 
##               85               50               50               50 
##     orangepurple           purple          redblue         redgreen 
##               50                9               49               51 
##          redpink        redpurple           silver            white 
##               56               55               26               48 
##        whiteblue        whitegold       whitegreen      whiteorange 
##               56               48               39               54 
##        whitepink      whitepurple         whitered      whiteyellow 
##               64               55               59               55 
##       yellowblue      yellowgreen     yelloworange       yellowpink 
##               54               51               54               58 
##     yellowpurple        yellowred 
##               54               57
hist(xtabs(~sl$beeCol[sl$trialNum == 1]))

mean(xtabs(~sl$beeCol[sl$trialNum == 1]))
## [1] 46.92857
sd(xtabs(~sl$beeCol[sl$trialNum == 1]))
## [1] 16.37725
# number of visits for second trial
xtabs(~sl$beeCol[sl$trialNum == 2])
## sl$beeCol[sl$trialNum == 2]
##             blue             gold          goldred            green 
##                0                0                0                4 
##             lime         limeblue         limegold        limegreen 
##                0               54                0                0 
##       limeorange         limepink       limepurple limepurpleyellow 
##               51                0              243                0 
##          limered       limesilver        limewhite       limeyellow 
##              281                0               51                0 
##           orange       orangeblue      orangegreen       orangepink 
##                0                0                0                0 
##     orangepurple           purple          redblue         redgreen 
##                0                0               60               53 
##          redpink        redpurple           silver            white 
##               58               90                0               67 
##        whiteblue        whitegold       whitegreen      whiteorange 
##              314               66                0              260 
##        whitepink      whitepurple         whitered      whiteyellow 
##              143                0              195               47 
##       yellowblue      yellowgreen     yelloworange       yellowpink 
##              128               50              176              142 
##     yellowpurple        yellowred 
##              346               56
# number of trials per bee and treatment
xtabs(~sl$beeCol +  sl$trt)
##                   sl$trt
## sl$beeCol          full high  low
##   blue               36    0    0
##   gold               54    0    0
##   goldred             3    0    0
##   green              13    0    0
##   lime               28    0    0
##   limeblue          530    0    0
##   limegold           54    0    0
##   limegreen          50    0    0
##   limeorange         84    0    0
##   limepink           51    0    0
##   limepurple        109    0  243
##   limepurpleyellow   58    0    0
##   limered            50    0 3424
##   limesilver          3    0    0
##   limewhite         101    0    0
##   limeyellow         52    0    0
##   orange             85  416    0
##   orangeblue         50    0    0
##   orangegreen        50    0    0
##   orangepink         50    0    0
##   orangepurple       50    0    0
##   purple              9    0    0
##   redblue           673    0    0
##   redgreen          530    0    0
##   redpink            56    0  621
##   redpurple          55  163    0
##   silver             26    0    0
##   white             283    0    0
##   whiteblue          56    0 2730
##   whitegold          48    0 1177
##   whitegreen         39    0    0
##   whiteorange        54    0 1068
##   whitepink         135 1049    0
##   whitepurple        55    0    0
##   whitered           59 1285    0
##   whiteyellow       157    0    0
##   yellowblue         54 1533    0
##   yellowgreen       532    0    0
##   yelloworange       54 2059    0
##   yellowpink         58    0 2251
##   yellowpurple       54    0 1189
##   yellowred         547    0    0
colSums(xtabs(~sl$beeCol +  sl$trt))
##  full  high   low 
##  5095  6505 12703
xtabs(~sl$beeCol +  sl$IT_imputed)
##                   sl$IT_imputed
## sl$beeCol          3.01  3.4 3.48 3.61 3.67 3.74 3.75 3.76 3.82 3.87 3.96
##   blue                0    0    0   36    0    0    0    0    0    0    0
##   gold                0    0    0    0    0    0    0    0    0    0    0
##   goldred             3    0    0    0    0    0    0    0    0    0    0
##   green               0    0    0    0    0    0    0    0    0    0    0
##   lime                0    0    0    0    0    0    0    0    0    0    0
##   limeblue            0    0    0    0    0    0    0    0    0    0    0
##   limegold            0    0    0    0    0    0    0    0    0    0    0
##   limegreen           0    0    0    0    0   50    0    0    0    0    0
##   limeorange          0    0    0    0    0    0    0    0    0    0    0
##   limepink            0    0    0    0    0    0    0    0    0    0    0
##   limepurple          0    0  352    0    0    0    0    0    0    0    0
##   limepurpleyellow    0    0    0    0    0    0   58    0    0    0    0
##   limered             0    0    0    0    0    0    0    0    0    0    0
##   limesilver          0    3    0    0    0    0    0    0    0    0    0
##   limewhite           0    0    0    0    0    0    0  101    0    0    0
##   limeyellow          0    0    0    0   52    0    0    0    0    0    0
##   orange              0    0    0    0    0    0    0    0    0    0    0
##   orangeblue          0    0    0    0    0    0    0    0    0    0    0
##   orangegreen         0    0    0    0    0    0    0    0    0    0    0
##   orangepink          0    0    0    0    0    0    0    0    0    0   50
##   orangepurple        0    0    0    0    0    0    0    0    0    0    0
##   purple              0    0    0    0    0    0    0    0    0    0    0
##   redblue             0    0    0    0    0    0  673    0    0    0    0
##   redgreen            0    0    0    0    0    0    0    0    0    0    0
##   redpink             0    0    0    0    0    0    0    0    0    0    0
##   redpurple           0    0    0    0    0    0    0    0    0    0    0
##   silver              0    0    0    0    0    0    0    0    0    0    0
##   white               0    0    0    0    0    0    0    0    0    0    0
##   whiteblue           0    0    0    0    0    0    0    0    0    0    0
##   whitegold           0    0    0    0    0    0    0    0    0    0    0
##   whitegreen          0    0    0    0    0    0    0    0    0    0    0
##   whiteorange         0    0    0    0    0    0    0    0    0    0    0
##   whitepink           0    0    0    0    0    0    0    0    0    0    0
##   whitepurple         0    0    0    0    0    0    0    0    0    0    0
##   whitered            0    0    0    0    0    0    0    0    0    0    0
##   whiteyellow         0    0    0    0    0    0    0    0    0    0    0
##   yellowblue          0    0    0    0    0    0    0    0    0    0    0
##   yellowgreen         0    0    0    0    0    0    0    0    0    0    0
##   yelloworange        0    0    0    0    0    0    0    0    0    0    0
##   yellowpink          0    0    0    0    0    0    0    0    0 2309    0
##   yellowpurple        0    0    0    0    0    0    0    0 1243    0    0
##   yellowred           0    0    0    0    0    0    0    0    0    0    0
##                   sl$IT_imputed
## sl$beeCol          3.98 3.99    4 4.02 4.03 4.05 4.07 4.12 4.13 4.14 4.15
##   blue                0    0    0    0    0    0    0    0    0    0    0
##   gold                0    0    0    0    0   54    0    0    0    0    0
##   goldred             0    0    0    0    0    0    0    0    0    0    0
##   green               0    0    0    0    0    0    0    0    0    0    0
##   lime                0    0    0   28    0    0    0    0    0    0    0
##   limeblue            0    0    0    0    0    0    0    0    0    0    0
##   limegold            0    0    0    0    0    0    0    0    0    0    0
##   limegreen           0    0    0    0    0    0    0    0    0    0    0
##   limeorange         84    0    0    0    0    0    0    0    0    0    0
##   limepink            0    0    0    0    0    0    0    0   51    0    0
##   limepurple          0    0    0    0    0    0    0    0    0    0    0
##   limepurpleyellow    0    0    0    0    0    0    0    0    0    0    0
##   limered             0    0    0    0    0    0    0    0    0    0 3474
##   limesilver          0    0    0    0    0    0    0    0    0    0    0
##   limewhite           0    0    0    0    0    0    0    0    0    0    0
##   limeyellow          0    0    0    0    0    0    0    0    0    0    0
##   orange              0  501    0    0    0    0    0    0    0    0    0
##   orangeblue          0    0    0    0    0    0    0    0    0    0    0
##   orangegreen         0    0    0    0    0    0    0    0    0    0    0
##   orangepink          0    0    0    0    0    0    0    0    0    0    0
##   orangepurple        0    0    0   50    0    0    0    0    0    0    0
##   purple              0    9    0    0    0    0    0    0    0    0    0
##   redblue             0    0    0    0    0    0    0    0    0    0    0
##   redgreen            0    0    0    0    0    0    0  530    0    0    0
##   redpink             0    0    0    0    0    0    0    0    0    0    0
##   redpurple           0    0    0    0    0    0    0    0    0    0    0
##   silver              0    0    0    0    0    0    0    0    0    0    0
##   white               0    0    0    0  283    0    0    0    0    0    0
##   whiteblue           0    0    0    0    0    0    0    0    0    0    0
##   whitegold           0    0    0    0    0    0    0    0    0    0    0
##   whitegreen          0    0   39    0    0    0    0    0    0    0    0
##   whiteorange         0    0    0    0    0    0    0    0    0    0    0
##   whitepink           0    0    0    0    0    0    0    0    0    0    0
##   whitepurple         0    0    0    0    0    0    0    0    0    0    0
##   whitered            0    0    0    0    0    0    0    0    0    0    0
##   whiteyellow         0    0    0    0    0    0    0    0    0    0    0
##   yellowblue          0    0    0    0    0    0    0    0    0    0    0
##   yellowgreen         0    0    0    0    0    0    0    0    0  532    0
##   yelloworange        0    0    0    0    0    0    0    0    0    0    0
##   yellowpink          0    0    0    0    0    0    0    0    0    0    0
##   yellowpurple        0    0    0    0    0    0    0    0    0    0    0
##   yellowred           0    0    0    0    0    0  547    0    0    0    0
##                   sl$IT_imputed
## sl$beeCol          4.21 4.23 4.24 4.31 4.34 4.37 4.42 4.45 4.46 4.59 4.61
##   blue                0    0    0    0    0    0    0    0    0    0    0
##   gold                0    0    0    0    0    0    0    0    0    0    0
##   goldred             0    0    0    0    0    0    0    0    0    0    0
##   green              13    0    0    0    0    0    0    0    0    0    0
##   lime                0    0    0    0    0    0    0    0    0    0    0
##   limeblue            0    0    0    0    0    0    0    0    0    0    0
##   limegold            0    0    0   54    0    0    0    0    0    0    0
##   limegreen           0    0    0    0    0    0    0    0    0    0    0
##   limeorange          0    0    0    0    0    0    0    0    0    0    0
##   limepink            0    0    0    0    0    0    0    0    0    0    0
##   limepurple          0    0    0    0    0    0    0    0    0    0    0
##   limepurpleyellow    0    0    0    0    0    0    0    0    0    0    0
##   limered             0    0    0    0    0    0    0    0    0    0    0
##   limesilver          0    0    0    0    0    0    0    0    0    0    0
##   limewhite           0    0    0    0    0    0    0    0    0    0    0
##   limeyellow          0    0    0    0    0    0    0    0    0    0    0
##   orange              0    0    0    0    0    0    0    0    0    0    0
##   orangeblue          0    0    0    0    0    0    0    0    0    0   50
##   orangegreen         0    0    0    0    0   50    0    0    0    0    0
##   orangepink          0    0    0    0    0    0    0    0    0    0    0
##   orangepurple        0    0    0    0    0    0    0    0    0    0    0
##   purple              0    0    0    0    0    0    0    0    0    0    0
##   redblue             0    0    0    0    0    0    0    0    0    0    0
##   redgreen            0    0    0    0    0    0    0    0    0    0    0
##   redpink           677    0    0    0    0    0    0    0    0    0    0
##   redpurple           0    0    0    0  218    0    0    0    0    0    0
##   silver              0   26    0    0    0    0    0    0    0    0    0
##   white               0    0    0    0    0    0    0    0    0    0    0
##   whiteblue           0    0    0    0 2786    0    0    0    0    0    0
##   whitegold           0    0    0    0    0    0    0    0    0    0    0
##   whitegreen          0    0    0    0    0    0    0    0    0    0    0
##   whiteorange         0    0    0    0    0    0    0    0 1122    0    0
##   whitepink           0    0    0    0    0    0    0    0    0 1184    0
##   whitepurple         0    0    0    0    0    0    0   55    0    0    0
##   whitered            0    0    0    0    0    0 1344    0    0    0    0
##   whiteyellow         0    0  157    0    0    0    0    0    0    0    0
##   yellowblue          0    0    0    0    0    0    0    0    0    0    0
##   yellowgreen         0    0    0    0    0    0    0    0    0    0    0
##   yelloworange        0    0    0    0 2113    0    0    0    0    0    0
##   yellowpink          0    0    0    0    0    0    0    0    0    0    0
##   yellowpurple        0    0    0    0    0    0    0    0    0    0    0
##   yellowred           0    0    0    0    0    0    0    0    0    0    0
##                   sl$IT_imputed
## sl$beeCol          4.63 4.69 4.94
##   blue                0    0    0
##   gold                0    0    0
##   goldred             0    0    0
##   green               0    0    0
##   lime                0    0    0
##   limeblue            0  530    0
##   limegold            0    0    0
##   limegreen           0    0    0
##   limeorange          0    0    0
##   limepink            0    0    0
##   limepurple          0    0    0
##   limepurpleyellow    0    0    0
##   limered             0    0    0
##   limesilver          0    0    0
##   limewhite           0    0    0
##   limeyellow          0    0    0
##   orange              0    0    0
##   orangeblue          0    0    0
##   orangegreen         0    0    0
##   orangepink          0    0    0
##   orangepurple        0    0    0
##   purple              0    0    0
##   redblue             0    0    0
##   redgreen            0    0    0
##   redpink             0    0    0
##   redpurple           0    0    0
##   silver              0    0    0
##   white               0    0    0
##   whiteblue           0    0    0
##   whitegold        1225    0    0
##   whitegreen          0    0    0
##   whiteorange         0    0    0
##   whitepink           0    0    0
##   whitepurple         0    0    0
##   whitered            0    0    0
##   whiteyellow         0    0    0
##   yellowblue          0    0 1587
##   yellowgreen         0    0    0
##   yelloworange        0    0    0
##   yellowpink          0    0    0
##   yellowpurple        0    0    0
##   yellowred           0    0    0
# hist(unique(sl$IT_imputed[sl$trt2 == "low"]))
# rug(unique(sl$IT_imputed[sl$trt2 == "low"]))
# 
# hist(unique(sl$IT_imputed[sl$trt2 == "high"]))
# rug(unique(sl$IT_imputed[sl$trt2 == "high"]))

# number of trials per bee and treatment
xtabs(~sl$beeCol +  sl$trialNum)
##                   sl$trialNum
## sl$beeCol            1   2   3   4   5   6   7   8   9  10  11  12
##   blue              36   0   0   0   0   0   0   0   0   0   0   0
##   gold              54   0   0   0   0   0   0   0   0   0   0   0
##   goldred            3   0   0   0   0   0   0   0   0   0   0   0
##   green              9   4   0   0   0   0   0   0   0   0   0   0
##   lime              28   0   0   0   0   0   0   0   0   0   0   0
##   limeblue          50  54  50  50  53  55  52  62  53  51   0   0
##   limegold          54   0   0   0   0   0   0   0   0   0   0   0
##   limegreen         50   0   0   0   0   0   0   0   0   0   0   0
##   limeorange        33  51   0   0   0   0   0   0   0   0   0   0
##   limepink          51   0   0   0   0   0   0   0   0   0   0   0
##   limepurple        53 243  56   0   0   0   0   0   0   0   0   0
##   limepurpleyellow  58   0   0   0   0   0   0   0   0   0   0   0
##   limered           50 281 351 298 486 537 393 170  96 105 707   0
##   limesilver         3   0   0   0   0   0   0   0   0   0   0   0
##   limewhite         50  51   0   0   0   0   0   0   0   0   0   0
##   limeyellow        52   0   0   0   0   0   0   0   0   0   0   0
##   orange            85   0  11  88  14  85 137  81   0   0   0   0
##   orangeblue        50   0   0   0   0   0   0   0   0   0   0   0
##   orangegreen       50   0   0   0   0   0   0   0   0   0   0   0
##   orangepink        50   0   0   0   0   0   0   0   0   0   0   0
##   orangepurple      50   0   0   0   0   0   0   0   0   0   0   0
##   purple             9   0   0   0   0   0   0   0   0   0   0   0
##   redblue           49  60  67 131  72  50  60  55  33  49  47   0
##   redgreen          51  53  51  53  51  60  54  54  52  51   0   0
##   redpink           56  58  71  57  62  62  76 119  62  54   0   0
##   redpurple         55  90  73   0   0   0   0   0   0   0   0   0
##   silver            26   0   0   0   0   0   0   0   0   0   0   0
##   white             48  67  46  40  57  25   0   0   0   0   0   0
##   whiteblue         56 314 307 190 539 537 105 103 276 111 248   0
##   whitegold         48  66  76 174 179 198 216 135  64  69   0   0
##   whitegreen        39   0   0   0   0   0   0   0   0   0   0   0
##   whiteorange       54 260 432 376   0   0   0   0   0   0   0   0
##   whitepink         64 143  71  99  78 165 122 127 105 118  92   0
##   whitepurple       55   0   0   0   0   0   0   0   0   0   0   0
##   whitered          59 195 169 278 109  89 112  87 138 108   0   0
##   whiteyellow       55  47  55   0   0   0   0   0   0   0   0   0
##   yellowblue        54 128  92 114 149 116 251 198 167 118 200   0
##   yellowgreen       51  50  56  50  58  52  56  51  54  54   0   0
##   yelloworange      54 176 205 475 113 256 111 141 252 110 220   0
##   yellowpink        58 142 103 140 235 187 278 111 366 266 372  51
##   yellowpurple      54 346 360 330 153   0   0   0   0   0   0   0
##   yellowred         57  56  58  55  51  57  55  56  50  52   0   0
# tabulate num bees in each trt
tdf <- data.frame(xtabs(~sl$beeCol +  sl$trt))
tdf2 <- acast(tdf, sl.beeCol~sl.trt, value.var="Freq")
tdf3 <- as.data.frame(tdf2 > 0)

trtdf<-table(interaction(tdf3$full, tdf3$high, tdf3$low))
names(trtdf) <- c("FULL ONLY", "FULL-HIGH", "FULL-LOW", "NA")
trtdf
## FULL ONLY FULL-HIGH  FULL-LOW        NA 
##        28         6         8         0
# get number of buzzes in each treatment
colSums(tdf2)
##  full  high   low 
##  5095  6505 12703
# get the number of bees for each treatment (every bee has "full" treatment)
apply(tdf2, MARGIN = 2, function(x) length(x[x>0]))
## full high  low 
##   42    6    8
head(sl)
##   beeCol hive meanFreq IT_imputed index freq     amp
## 1   blue    4 395.2778       3.61     1  400 0.06328
## 2   blue    4 395.2778       3.61     2  360 0.34299
## 3   blue    4 395.2778       3.61     3  430 0.44099
## 4   blue    4 395.2778       3.61     4  420 0.16841
## 5   blue    4 395.2778       3.61     5  420 0.15284
## 6   blue    4 395.2778       3.61     6  410 0.18302
##                    datetime rewNum rewTF lowRewAmp highrewAmp
## 1  2016_11_22__08_23_47_721      1     T         0          5
## 2  2016_11_22__08_23_49_677      2     T         0          5
## 3  2016_11_22__08_23_50_139      3     T         0          5
## 4  2016_11_22__08_24_10_938      4     T         0          5
## 5  2016_11_22__08_24_18_013      5     T         0          5
## 6  2016_11_22__08_24_18_464      6     T         0          5
##                         BeeNumCol
## 1 BeeBlue_22Nov2016_Hive4_initial
## 2 BeeBlue_22Nov2016_Hive4_initial
## 3 BeeBlue_22Nov2016_Hive4_initial
## 4 BeeBlue_22Nov2016_Hive4_initial
## 5 BeeBlue_22Nov2016_Hive4_initial
## 6 BeeBlue_22Nov2016_Hive4_initial
##                                     accFile trialNum
## 1 2016_11_22__08_23_47_721_220_450_test.txt        1
## 2 2016_11_22__08_23_49_677_220_450_test.txt        1
## 3 2016_11_22__08_23_50_139_220_450_test.txt        1
## 4 2016_11_22__08_24_10_938_220_450_test.txt        1
## 5 2016_11_22__08_24_18_013_220_450_test.txt        1
## 6 2016_11_22__08_24_18_464_220_450_test.txt        1
##                 datetime_str lowFrq highFrq   IT  trt   amp_acc
## 1 2016-11-22 08:23:47.721000    220     450 3.61 full  6.222222
## 2 2016-11-22 08:23:49.677000    220     450 3.61 full 33.725664
## 3 2016-11-22 08:23:50.139000    220     450 3.61 full 43.361849
## 4 2016-11-22 08:24:10.938000    220     450 3.61 full 16.559489
## 5 2016-11-22 08:24:18.013000    220     450 3.61 full 15.028515
## 6 2016-11-22 08:24:18.464000    220     450 3.61 full 17.996067
dim(sl) # should be 24303 rows
## [1] 24303    21
# hive 5 is most common
table(sl$hive, useNA = 'always')
## 
##     3     4     5  <NA> 
##   513  1047 22743     0
## make sure all bee colors are lowercase
sl$beeCol <- tolower(sl$beeCol)

# make sure there are values lower than 220 and higher than 450 
# (the cutoff for buzzes used in the experiment)
hist(sl$freq, breaks = seq(215, 450, by = 10))

nrow(sl[sl$freq < 220 | sl$freq > 450,]) # should have 0 rows
## [1] 0
# look at treatments
xtabs(~sl$beeCol+ trt, data = sl )
##                   trt
## sl$beeCol          full high  low
##   blue               36    0    0
##   gold               54    0    0
##   goldred             3    0    0
##   green              13    0    0
##   lime               28    0    0
##   limeblue          530    0    0
##   limegold           54    0    0
##   limegreen          50    0    0
##   limeorange         84    0    0
##   limepink           51    0    0
##   limepurple        109    0  243
##   limepurpleyellow   58    0    0
##   limered            50    0 3424
##   limesilver          3    0    0
##   limewhite         101    0    0
##   limeyellow         52    0    0
##   orange             85  416    0
##   orangeblue         50    0    0
##   orangegreen        50    0    0
##   orangepink         50    0    0
##   orangepurple       50    0    0
##   purple              9    0    0
##   redblue           673    0    0
##   redgreen          530    0    0
##   redpink            56    0  621
##   redpurple          55  163    0
##   silver             26    0    0
##   white             283    0    0
##   whiteblue          56    0 2730
##   whitegold          48    0 1177
##   whitegreen         39    0    0
##   whiteorange        54    0 1068
##   whitepink         135 1049    0
##   whitepurple        55    0    0
##   whitered           59 1285    0
##   whiteyellow       157    0    0
##   yellowblue         54 1533    0
##   yellowgreen       532    0    0
##   yelloworange       54 2059    0
##   yellowpink         58    0 2251
##   yellowpurple       54    0 1189
##   yellowred         547    0    0
# find percentage reward by treatment
mean(grepl("[tT]", as.character(sl$rewTF))) # overall mean
## [1] 0.4155454
# percentage that were rewarded by treatment
tapply((grepl("[tT]", as.character(sl$rewTF))), INDEX= sl$trt, mean)
##      full      high       low 
## 1.0000000 0.3535742 0.2128631
# total number of trials for each treatment
tapply((grepl("[tT]", as.character(sl$rewTF))), INDEX= sl$trt, length)
##  full  high   low 
##  5095  6505 12703
# total number of rewards per treatment
tapply((grepl("[tT]", as.character(sl$rewTF))), INDEX= sl$trt, FUN = function(x) sum(x))
## full high  low 
## 5095 2300 2704
# total number of trials that were unrewarded per treatment
tapply((grepl("[tT]", as.character(sl$rewTF))), INDEX= sl$trt, FUN = function(x) sum(!x))
## full high  low 
##    0 4205 9999
# calculate mean freq per bee per trial
d1 <- as.data.frame((tapply(sl$freq, INDEX = list(sl$beeCol, sl$trt), mean)))
d2 <- as.data.frame((tapply(sl$freq, INDEX = list(sl$beeCol, sl$trt), length)))
d3 <- data.frame(d1, d2)
d3$beeCol = row.names(d3)


d4 <- melt(d3, id.vars=c("beeCol"), measure.vars=c("full", "high", "low" ), value.name = "buzzFreq")
d5 <- melt(d3, id.vars=c("beeCol"), measure.vars=c("full.1",  "high.1", "low.1" ), value.name = "NumBuzzes")


d6 <- data.frame(d4, d5)
d7 <- d6[!is.na(d6$buzzFreq),]
d7
##               beeCol variable buzzFreq         beeCol.1 variable.1
## 1               blue     full 395.2778             blue     full.1
## 2               gold     full 348.3333             gold     full.1
## 3            goldred     full 406.6667          goldred     full.1
## 4              green     full 308.4615            green     full.1
## 5               lime     full 374.6429             lime     full.1
## 6           limeblue     full 304.1509         limeblue     full.1
## 7           limegold     full 344.6296         limegold     full.1
## 8          limegreen     full 370.8000        limegreen     full.1
## 9         limeorange     full 321.4286       limeorange     full.1
## 10          limepink     full 308.4314         limepink     full.1
## 11        limepurple     full 340.6422       limepurple     full.1
## 12  limepurpleyellow     full 336.7241 limepurpleyellow     full.1
## 13           limered     full 354.6000          limered     full.1
## 14        limesilver     full 306.6667       limesilver     full.1
## 15         limewhite     full 309.7030        limewhite     full.1
## 16        limeyellow     full 339.4231       limeyellow     full.1
## 17            orange     full 362.4706           orange     full.1
## 18        orangeblue     full 301.4000       orangeblue     full.1
## 19       orangegreen     full 296.2000      orangegreen     full.1
## 20        orangepink     full 286.6000       orangepink     full.1
## 21      orangepurple     full 329.0000     orangepurple     full.1
## 22            purple     full 382.2222           purple     full.1
## 23           redblue     full 341.2927          redblue     full.1
## 24          redgreen     full 326.9057         redgreen     full.1
## 25           redpink     full 288.3929          redpink     full.1
## 26         redpurple     full 311.4545        redpurple     full.1
## 27            silver     full 305.3846           silver     full.1
## 28             white     full 319.1873            white     full.1
## 29         whiteblue     full 318.2143        whiteblue     full.1
## 30         whitegold     full 315.2083        whitegold     full.1
## 31        whitegreen     full 300.7692       whitegreen     full.1
## 32       whiteorange     full 324.4444      whiteorange     full.1
## 33         whitepink     full 333.0370        whitepink     full.1
## 34       whitepurple     full 328.1818      whitepurple     full.1
## 35          whitered     full 288.8136         whitered     full.1
## 36       whiteyellow     full 316.8153      whiteyellow     full.1
## 37        yellowblue     full 313.7037       yellowblue     full.1
## 38       yellowgreen     full 284.7368      yellowgreen     full.1
## 39      yelloworange     full 268.1481     yelloworange     full.1
## 40        yellowpink     full 323.6207       yellowpink     full.1
## 41      yellowpurple     full 354.4444     yellowpurple     full.1
## 42         yellowred     full 337.4589        yellowred     full.1
## 59            orange     high 381.5385           orange     high.1
## 68         redpurple     high 336.4417        redpurple     high.1
## 75         whitepink     high 326.2917        whitepink     high.1
## 77          whitered     high 321.5019         whitered     high.1
## 79        yellowblue     high 309.4455       yellowblue     high.1
## 81      yelloworange     high 311.5590     yelloworange     high.1
## 95        limepurple      low 344.2798       limepurple      low.1
## 97           limered      low 364.7138          limered      low.1
## 109          redpink      low 290.2576          redpink      low.1
## 113        whiteblue      low 348.4505        whiteblue      low.1
## 114        whitegold      low 329.6092        whitegold      low.1
## 116      whiteorange      low 361.1610      whiteorange      low.1
## 124       yellowpink      low 346.0640       yellowpink      low.1
## 125     yellowpurple      low 361.1102     yellowpurple      low.1
##     NumBuzzes
## 1          36
## 2          54
## 3           3
## 4          13
## 5          28
## 6         530
## 7          54
## 8          50
## 9          84
## 10         51
## 11        109
## 12         58
## 13         50
## 14          3
## 15        101
## 16         52
## 17         85
## 18         50
## 19         50
## 20         50
## 21         50
## 22          9
## 23        673
## 24        530
## 25         56
## 26         55
## 27         26
## 28        283
## 29         56
## 30         48
## 31         39
## 32         54
## 33        135
## 34         55
## 35         59
## 36        157
## 37         54
## 38        532
## 39         54
## 40         58
## 41         54
## 42        547
## 59        416
## 68        163
## 75       1049
## 77       1285
## 79       1533
## 81       2059
## 95        243
## 97       3424
## 109       621
## 113      2730
## 114      1177
## 116      1068
## 124      2251
## 125      1189
# this plot shows how individual bees changed
ggplot(d7, aes(x = variable, y = buzzFreq)) + 
  geom_boxplot(aes(fill = variable == "full"), alpha = 0.2) + 
  geom_point(aes(color = variable == "full", size = NumBuzzes)) + 
  geom_line(aes(group = beeCol), alpha = 0.2) + 
  scale_fill_viridis(name = "hihi", guide = FALSE, discrete = TRUE, option = "A", end = 0.6) + 
  scale_color_viridis(name = "hihi", guide = FALSE, discrete = TRUE, option = "A", end = 0.6) + 
  labs(x = "treatment")

Visualizations and modeling

# calculate trial averages and plot
sl$colNum = paste(sl$beeCol, sl$hive, sep = "_")
sl$trt <- as.character(sl$trt)
sl$trt[sl$trt == "full" & sl$trialNum >1 ] <- "full_2"


aggdata <- aggregate(sl$freq, by=list(colNum = sl$colNum,trialNum =sl$trialNum, trt = sl$trt), FUN=mean, na.rm=TRUE)
colnames(aggdata)[colnames(aggdata) == "x"] = "freq"

aggdata_sd <- aggregate(sl$freq, by=list(colNum = sl$colNum,trialNum =sl$trialNum, trt = sl$trt), FUN=sd, na.rm=TRUE)
colnames(aggdata_sd)[colnames(aggdata_sd) == "x"] = "freq_sd"

aggdata = merge(aggdata, aggdata_sd)
#aggdata$trt = sapply(1:nrow(aggdata), FUN = function(x){sl[sl$colNum == aggdata$colNum[x] & sl$trialNum == aggdata$trialNum[x], "trt"][1]})
aggdata = aggdata[order(aggdata$colNum, aggdata$trialNum, decreasing = FALSE), ]
agg_sm = aggdata[aggdata$trialNum <= 2, ]
rownames(agg_sm) = 1:nrow(agg_sm)
agg_sm
##                colNum trialNum    trt     freq  freq_sd
## 1              blue_4        1   full 395.2778 39.38717
## 2              gold_3        1   full 348.3333 37.15089
## 3           goldred_4        1   full 406.6667 32.14550
## 4             green_4        1   full 318.8889 57.32461
## 5             green_4        2 full_2 285.0000 59.16080
## 6              lime_5        1   full 374.6429 36.36055
## 7          limeblue_5        1   full 302.4000 34.73280
## 8          limeblue_5        2 full_2 299.8148 39.02059
## 9          limegold_5        1   full 344.6296 44.07168
## 10        limegreen_5        1   full 370.8000 29.40498
## 11       limeorange_5        1   full 285.7576 35.26953
## 12       limeorange_5        2 full_2 344.5098 26.70683
## 13         limepink_5        1   full 308.4314 43.23760
## 14       limepurple_5        1   full 352.2642 38.66232
## 15       limepurple_5        2    low 344.2798 32.55948
## 16 limepurpleyellow_5        1   full 336.7241 50.76219
## 17          limered_5        1   full 354.6000 37.91559
## 18          limered_5        2    low 365.3381 36.20344
## 19       limesilver_5        1   full 306.6667 11.54701
## 20        limewhite_5        1   full 292.0000 44.90352
## 21        limewhite_5        2 full_2 327.0588 35.51305
## 22       limeyellow_5        1   full 339.4231 39.37818
## 23           orange_3        1   full 388.5294 38.22837
## 24           orange_5        1   full 345.0980 33.60789
## 25       orangeblue_5        1   full 301.4000 45.17585
## 26      orangegreen_5        1   full 296.2000 40.45103
## 27       orangepink_5        1   full 286.6000 27.15113
## 28     orangepurple_5        1   full 329.0000 25.65469
## 29           purple_3        1   full 382.2222 21.08185
## 30          redblue_4        1   full 342.2449 40.53113
## 31          redblue_4        2 full_2 335.6667 51.89189
## 32         redgreen_5        1   full 334.7059 34.19666
## 33         redgreen_5        2 full_2 314.5283 31.53545
## 34          redpink_5        1   full 288.3929 41.06654
## 35          redpink_5        2    low 275.3448 34.24450
## 36        redpurple_5        1   full 311.4545 31.64726
## 37        redpurple_5        2   high 336.4444 40.75995
## 38           silver_5        1   full 305.3846 38.39070
## 39            white_4        1   full 343.1250 47.49720
## 40            white_4        2 full_2 342.2388 38.95689
## 41        whiteblue_5        1   full 318.2143 27.17667
## 42        whiteblue_5        2    low 351.6242 29.13402
## 43        whitegold_5        1   full 315.2083 34.20772
## 44        whitegold_5        2    low 307.2727 45.76009
## 45       whitegreen_4        1   full 300.7692 33.43372
## 46      whiteorange_5        1   full 324.4444 51.71280
## 47      whiteorange_5        2    low 355.2308 44.74234
## 48        whitepink_5        1   full 336.8750 42.68285
## 49        whitepink_5        2   high 323.9161 28.55770
## 50      whitepurple_5        1   full 328.1818 43.50781
## 51         whitered_5        1   full 288.8136 36.05956
## 52         whitered_5        2   high 304.6667 42.10035
## 53      whiteyellow_5        1   full 303.4545 39.07344
## 54      whiteyellow_5        2 full_2 312.5532 36.56252
## 55       yellowblue_5        1   full 313.7037 44.86021
## 56       yellowblue_5        2   high 327.9688 49.76099
## 57      yellowgreen_5        1   full 298.0392 40.64577
## 58      yellowgreen_5        2 full_2 300.6000 27.58364
## 59     yelloworange_5        1   full 268.1481 36.18996
## 60     yelloworange_5        2   high 313.9205 39.62646
## 61       yellowpink_5        1   full 323.6207 39.98676
## 62       yellowpink_5        2    low 335.7042 41.00507
## 63     yellowpurple_5        1   full 354.4444 36.37678
## 64     yellowpurple_5        2    low 373.0058 34.17237
## 65        yellowred_5        1   full 320.1754 41.85396
## 66        yellowred_5        2 full_2 333.0357 29.96481
agg_sm[duplicated(data.frame(agg_sm$colNum, agg_sm$trialNum)), ]
## [1] colNum   trialNum trt      freq     freq_sd 
## <0 rows> (or 0-length row.names)
diffdf <- sapply(unique(agg_sm$colNum), FUN = function(x){
  tmp = agg_sm[agg_sm$colNum == x, ]
  if(nrow(tmp) <= 1)
    diff = NA
  else
    diff = tmp$freq[tmp$trialNum == 2] - tmp$freq[tmp$trialNum == 1]
  return(diff)
})

trtDF = sapply(unique(agg_sm$colNum), FUN = function(x){
  tmp = agg_sm[agg_sm$colNum == x, ]
  ttrs = paste(tmp$trt[tmp$trialNum == 1], tmp$trt[tmp$trialNum == 2], sep = "_")
  return(ttrs)
})

buzzdiffs = data.frame(trtDF, diffdf)

ggplot(buzzdiffs, aes(x = trtDF, y= diffdf)) + 
  geom_boxplot() + 
  geom_point()
## Warning: Removed 20 rows containing non-finite values (stat_boxplot).
## Warning: Removed 20 rows containing missing values (geom_point).

agg2 = aggregate(sl$freq, by=list(colNum = sl$colNum, fullTrt = sl$trt == "full", trt = sl$trt), FUN=mean, na.rm=TRUE)
colnames(agg2)[colnames(agg2) == "x"] = "freq"
agg2$trt = as.character(agg2$trt)

diffdf <- t(as.data.frame(t(sapply(unique(agg2$colNum), FUN = function(x){
  tmp = agg2[agg2$colNum == x, ]
  if(nrow(tmp) <= 1)
    return(NA)
  if (length(unique(tmp$trt)) > 2){
    tmp = tmp[tmp$trt != "full_2", ]
  }
  diff =  tmp$freq[!tmp$fullTrt] - tmp$freq[tmp$fullTrt]
  return(diff)

}))))

length(diffdf)
## [1] 43
trtDF = sapply(unique(agg2$colNum), FUN = function(x){
  tmp = agg2[agg2$colNum == x, ]
  if (length(unique(tmp$trt)) > 2){
    tmp = tmp[tmp$trt != "full_2", ]
  }
  ttrs = paste(tmp$trt[tmp$fullTrt], tmp$trt[!tmp$fullTrt], sep = "_")
  return(ttrs)
})

length(trtDF)
## [1] 43
buzzdiffs = data.frame(trtDF, diffdf)
tapply(buzzdiffs$diffdf, INDEX = buzzdiffs$trtDF, mean)
##       full_ full_full_2   full_high    full_low 
##          NA    4.851689   13.209001   14.307123
ggplot(droplevels(buzzdiffs[buzzdiffs$trtDF != "full_", ]), aes(x = trtDF, y= as.numeric(diffdf))) + 
  geom_boxplot() + 
  geom_point() + 
  labs(y = "diff in frequency -- averages")

MODELING

sl$trt = relevel(factor(sl$trt), ref = "full")


sl$trt2 = mapvalues(sl$trt, from = c("full", "full_2", "high", "low"), 
                    to = c("full", "full", "high", "low"))
# summary for paper
# fit a varying slope and intercept for colNum (bee ID), and allow the slope of the trialNum
# variable to vary by colNum (beeID)
m2 = lmer(freq ~ trt2 + hive + trialNum + IT_imputed + (1+trialNum|colNum), data = sl, REML = FALSE)

# this model estimates a global intercept
# random effect intercept for colNum (beeID)
# a single global estimate for trialNum
# the effect of trialNum within each level of colNum
# the correlation between intercept of trialNum across levels of colNum


vif.mer <- function (fit) {
    ## adapted from rms::vif

    v <- vcov(fit)
    nam <- names(fixef(fit))

    ## exclude intercepts
    ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)"))
    if (ns > 0) {
        v <- v[-(1:ns), -(1:ns), drop = FALSE]
        nam <- nam[-(1:ns)]
    }

    d <- diag(v)^0.5
    v <- diag(solve(v/(d %o% d)))
    names(v) <- nam
    v
}

vif.mer(m2) # fine
##   trt2high    trt2low      hive4      hive5   trialNum IT_imputed 
##   1.010841   1.004759   2.648702   2.664859   1.014513   1.143183
m2_1 = lmer(freq ~ trt2* IT_imputed + hive + trialNum  + (1+trialNum|colNum), data = sl, REML = FALSE)

# compare models with BIC
BIC(m2, m2_1) # m2 (no interaction) is better
##      df      BIC
## m2   11 246086.5
## m2_1 13 246090.2
anova(m2, m2_1) # Note: this disagrees with BIC
## Data: sl
## Models:
## m2: freq ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum | 
## m2:     colNum)
## m2_1: freq ~ trt2 * IT_imputed + hive + trialNum + (1 + trialNum | 
## m2_1:     colNum)
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m2   11 245997 246087 -122988   245975                             
## m2_1 13 245985 246090 -122979   245959 16.547      2  0.0002552 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3 = update(m2, .~. - hive)
BIC(m2, m3) # m3 is better (without hive)
##    df      BIC
## m2 11 246086.5
## m3  9 246073.9
anova(m2, m3) # Note: disagrees with BIC
## Data: sl
## Models:
## m3: freq ~ trt2 + trialNum + IT_imputed + (1 + trialNum | colNum)
## m2: freq ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum | 
## m2:     colNum)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m3  9 246001 246074 -122991   245983                           
## m2 11 245997 246087 -122988   245975 7.5417      2    0.02303 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# refit best model with REML = TRUE
m3 <- update(m3, .~., REML = TRUE)
summary(m3) # summary for paper -- leave trial number in, to account for random effect
## Linear mixed model fit by REML ['lmerMod']
## Formula: freq ~ trt2 + trialNum + IT_imputed + (1 + trialNum | colNum)
##    Data: sl
## 
## REML criterion at convergence: 245962
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0242 -0.5161  0.1663  0.6809  3.9135 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  colNum   (Intercept) 1072     32.74         
##           trialNum     112     10.58    -0.73
##  Residual             1439     37.93         
## Number of obs: 24303, groups:  colNum, 43
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  479.138     44.361  10.801
## trt2high      11.161      2.294   4.865
## trt2low       17.718      1.916   9.248
## trialNum       1.448      2.133   0.679
## IT_imputed   -37.291     10.773  -3.462
## 
## Correlation of Fixed Effects:
##            (Intr) trt2hg trt2lw trilNm
## trt2high    0.069                     
## trt2low    -0.005  0.009              
## trialNum    0.007 -0.042 -0.037       
## IT_imputed -0.993 -0.072  0.000 -0.093
plot(m3)

qqnorm(ranef(m3)$colNum[[1]])
qqline(ranef(m3)$colNum[[1]])

sjp.lmer(m3) # plot random effects to find any outliers
## Plotting random effects...

sjp.lmer(m3, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects
## Computing p-values via Wald-statistics approximation (treating t as Wald z).

# post-hoc tests -- pvals for paper
summary(glht(m3, linfct = mcp(trt2 = "Tukey")), test = adjusted("none"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = freq ~ trt2 + trialNum + IT_imputed + (1 + trialNum | 
##     colNum), data = sl, REML = TRUE)
## 
## Linear Hypotheses:
##                  Estimate Std. Error z value Pr(>|z|)    
## high - full == 0   11.161      2.294   4.865 1.14e-06 ***
## low - full == 0    17.718      1.916   9.248  < 2e-16 ***
## low - high == 0     6.557      2.975   2.204   0.0275 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)
## bonf adjusted pvals
summary(glht(m3, linfct = mcp(trt2 = "Tukey")), test = adjusted("bonf"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = freq ~ trt2 + trialNum + IT_imputed + (1 + trialNum | 
##     colNum), data = sl, REML = TRUE)
## 
## Linear Hypotheses:
##                  Estimate Std. Error z value Pr(>|z|)    
## high - full == 0   11.161      2.294   4.865 3.43e-06 ***
## low - full == 0    17.718      1.916   9.248  < 2e-16 ***
## low - high == 0     6.557      2.975   2.204   0.0826 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)

Generate CI’s

# set number of bOotstrap samples
nbootSims = 1000

ITmean = mean(tapply(sl$IT_imputed, INDEX = sl$beeCol, FUN = function(x) x[1] ))
print(ITmean)
## [1] 4.092619
# don't need hive, because that's not in the model we chose (above)
pframe <- data.frame(trt2 = levels(droplevels(sl$trt2)), IT_imputed = ITmean,  colNum = 99999, trialNum = 2)
pframe$freq <- 0
pp <- predict(m3, newdata = pframe, re.form=NA, type = 'response') # re.form sets all random effects to 0


### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(m3, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims)
print(paste("Number of bootstrap samples", nrow(bb2$t)))
## [1] "Number of bootstrap samples 1000"
bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
pframe$blo<-bb2_se[1,]
pframe$bhi<-bb2_se[2,]
pframe$predMean <- pp
pframe <- pframe[, c('trt2', "blo", "bhi", "predMean")]
pframe
##   trt2      blo      bhi predMean
## 1 full 321.7666 336.2141 329.4150
## 2 high 332.6790 348.3194 340.5761
## 3  low 339.1716 355.1462 347.1330

Make plots for paper

#plot 95% confidence intervals
# "Mean and bootstrap CI based on fixed-effects uncertainty ONLY"
# rename labels for plot
pframe$trt3 = mapvalues(pframe$trt2, from = c("full", "high", "low"), 
                        to = c("Full range\n(220 - 450 Hz)", "High range\n(340 - 390 Hz)", "Low range\n(220 - 330 Hz)"))

pframe
##   trt2      blo      bhi predMean                       trt3
## 1 full 321.7666 336.2141 329.4150 Full range\n(220 - 450 Hz)
## 2 high 332.6790 348.3194 340.5761 High range\n(340 - 390 Hz)
## 3  low 339.1716 355.1462 347.1330  Low range\n(220 - 330 Hz)
g0 <- ggplot(pframe, aes(x=trt3, y=predMean))+
     geom_point()+ 
     labs(y = "Sonication frequency (Hz)", x = "Frequency range for reward") + 
     geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1)+
  theme(plot.background = element_rect(fill = "transparent",colour = NA), 
        panel.background = element_rect(fill = "transparent",colour = NA)) + 
  annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 360, label=c("a", "b", "c"),
                color="black") 
g0

ggsave(plot = g0, filename = file.path(figDir, "SonicationFreqPredsAndCI_unadjusted.png"), width = 5, height = 4)

# save with transparent background
g0 <- ggplot(pframe, aes(x=trt3, y=predMean))+
     geom_point()+ 
     labs(y = "Sonication frequency (Hz)", x = "Frequency range for reward") + 
     geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1)+
  theme(plot.background = element_rect(fill = "transparent",colour = NA), 
        panel.background = element_rect(fill = "transparent",colour = NA)) + 
  annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 360, label=c("a", "b", "b"),
                color="black") 
g0

ggsave(plot = g0, filename = file.path(figDir, "SonicationFreqPredsAndCI_adjustedPvals.pdf"), width = 5, height = 4)




# use the "effects" library to get CI's
# this basically gives the same result as bootstrapping
# it's a good double check
ee <- as.data.frame(Effect(c("trt2"), 
                           fixed.predictors =list(given.values = c(trialNum = 2, 
                                                                   IT_imputed = ITmean)),  
                           m3) )


# compare two methods -- very similar
ee
##   trt2      fit       se    lower    upper
## 1 full 329.4150 3.798940 321.9689 336.8612
## 2 high 340.5761 4.280727 332.1856 348.9666
## 3  low 347.1330 4.081495 339.1330 355.1330
pframe
##   trt2      blo      bhi predMean                       trt3
## 1 full 321.7666 336.2141 329.4150 Full range\n(220 - 450 Hz)
## 2 high 332.6790 348.3194 340.5761 High range\n(340 - 390 Hz)
## 3  low 339.1716 355.1462 347.1330  Low range\n(220 - 330 Hz)
g1 <- ggplot(ee, aes(x=trt2, y=fit))+
     geom_point()+ 
     labs(y = "Sonication Frequency (Hz)", x = "Treatment") + 
     geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.1)+
     theme_classic() + 
    annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 360, label=c("a", "b", "c"),
                color="black") 
g1

____ END OF ANALYSIS FOR PAPER ____

Analysis for amplitude

# refref: here try log amp_acc and REML = FALSE with BIC

# interaction model
# note log-transformation to make model fit assumptions better
maa0 = lmer(log(amp_acc) ~ trt2* IT_imputed + hive + trialNum  + (1+trialNum|colNum), data = sl, REML = FALSE)

# main effect model
maa1 = lmer(log(amp_acc) ~ trt2 + hive + trialNum + IT_imputed + (1+trialNum|colNum), data = sl, REML = FALSE)
BIC(maa0, maa1) # use no interaction (keep maa1)
##      df      BIC
## maa0 13 44377.70
## maa1 11 44360.43
anova(maa0, maa1) # note that this agrees with likelihood ratio test
## Data: sl
## Models:
## maa1: log(amp_acc) ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum | 
## maa1:     colNum)
## maa0: log(amp_acc) ~ trt2 * IT_imputed + hive + trialNum + (1 + trialNum | 
## maa0:     colNum)
##      Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## maa1 11 44271 44360 -22125    44249                         
## maa0 13 44272 44378 -22123    44246 2.9209      2     0.2321
maa2 = update(maa1, .~. - trt2)
BIC(maa1, maa2) # keep treatment (maa1)
##      df      BIC
## maa1 11 44360.43
## maa2  9 44375.90
anova(maa1, maa2) # agrees with LRT; p-val for paper, if needed
## Data: sl
## Models:
## maa2: log(amp_acc) ~ hive + trialNum + IT_imputed + (1 + trialNum | 
## maa2:     colNum)
## maa1: log(amp_acc) ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum | 
## maa1:     colNum)
##      Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
## maa2  9 44303 44376 -22142    44285                             
## maa1 11 44271 44360 -22125    44249 35.665      2  1.801e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
maa3 <- update(maa1, .~. - hive)
BIC(maa1, maa3) # remove hive (maa3)
##      df      BIC
## maa1 11 44360.43
## maa3  9 44345.85
anova(maa1, maa3) # agrees with LRT
## Data: sl
## Models:
## maa3: log(amp_acc) ~ trt2 + trialNum + IT_imputed + (1 + trialNum | 
## maa3:     colNum)
## maa1: log(amp_acc) ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum | 
## maa1:     colNum)
##      Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)  
## maa3  9 44273 44346 -22128    44255                           
## maa1 11 44271 44360 -22125    44249 5.6143      2    0.06038 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# refit with REML = TRUE for paper
maa3 <- update(maa3, .~., REML = TRUE)
summary(maa3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(amp_acc) ~ trt2 + trialNum + IT_imputed + (1 + trialNum |  
##     colNum)
##    Data: sl
## 
## REML criterion at convergence: 44279.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8377 -0.5961  0.1156  0.7367  2.9119 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  colNum   (Intercept) 0.07079  0.26607       
##           trialNum    0.00203  0.04506  -0.57
##  Residual             0.35872  0.59893       
## Number of obs: 24303, groups:  colNum, 43
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 2.398060   0.455225   5.268
## trt2high    0.194526   0.034446   5.647
## trt2low     0.067759   0.029326   2.311
## trialNum    0.005513   0.010107   0.545
## IT_imputed  0.324122   0.110268   2.939
## 
## Correlation of Fixed Effects:
##            (Intr) trt2hg trt2lw trilNm
## trt2high    0.107                     
## trt2low    -0.009  0.029              
## trialNum    0.036 -0.089 -0.083       
## IT_imputed -0.995 -0.114  0.000 -0.087
# diagnostics
plot(maa3)

qqnorm(ranef(maa1)$colNum[[1]])
qqline(ranef(maa1)$colNum[[1]])

sjp.lmer(maa3) # plot random effects to find any outliers
## Plotting random effects...

sjp.lmer(maa3, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects
## Computing p-values via Wald-statistics approximation (treating t as Wald z).

# post-hoc tests -- pvals for paper
summary(glht(maa3, linfct = mcp(trt2 = "Tukey")), test = adjusted("none"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = log(amp_acc) ~ trt2 + trialNum + IT_imputed + 
##     (1 + trialNum | colNum), data = sl, REML = TRUE)
## 
## Linear Hypotheses:
##                  Estimate Std. Error z value Pr(>|z|)    
## high - full == 0  0.19453    0.03445   5.647 1.63e-08 ***
## low - full == 0   0.06776    0.02933   2.311  0.02086 *  
## low - high == 0  -0.12677    0.04459  -2.843  0.00447 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)
# post-hoc tests with bonf adjustment
summary(glht(maa3, linfct = mcp(trt2 = "Tukey")), test = adjusted("bonf"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = log(amp_acc) ~ trt2 + trialNum + IT_imputed + 
##     (1 + trialNum | colNum), data = sl, REML = TRUE)
## 
## Linear Hypotheses:
##                  Estimate Std. Error z value Pr(>|z|)    
## high - full == 0  0.19453    0.03445   5.647 4.89e-08 ***
## low - full == 0   0.06776    0.02933   2.311   0.0626 .  
## low - high == 0  -0.12677    0.04459  -2.843   0.0134 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)

Generate CI’s for amplitude

# set number of bootstrap samples
nbootSims2 = 1000

# don't need to include hive, but will set trialNum to 2
pframe <- data.frame(trt2 = levels(droplevels(sl$trt2)), 
                     IT_imputed = ITmean,  colNum = 99999, trialNum = 2)
pframe$amp_acc <- 0
# exponentiate to get back to original scale
pp <- exp(predict(maa3, newdata = pframe, re.form=NA, type = 'response')) # re.form sets all random effects to 0


### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(maa3, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims2)
print(paste("Number of bootstrap samples", nrow(bb2$t)))
## [1] "Number of bootstrap samples 1000"
bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))

#exponentiate to get back to orignal scale
pframe$blo<-exp(bb2_se[1,])
pframe$bhi<-exp(bb2_se[2,])
pframe$predMean <- pp
pframe <- pframe[, c('trt2', "blo", "bhi", "predMean")]
pframe
##   trt2      blo      bhi predMean
## 1 full 39.05177 45.01982 41.91285
## 2 high 46.44376 56.04234 50.91299
## 3  low 41.17250 48.98075 44.85125

Make plots for amplitude for paper

#plot 95% confidence intervals
# "Mean and bootstrap CI based on fixed-effects uncertainty ONLY"
pframe$trt3 = mapvalues(pframe$trt2, from = c("full", "high", "low"), 
                        to = c("Full range\n(220 - 450 Hz)", "High range\n(340 - 390 Hz)", "Low range\n(220 - 330 Hz)"))

pframe
##   trt2      blo      bhi predMean                       trt3
## 1 full 39.05177 45.01982 41.91285 Full range\n(220 - 450 Hz)
## 2 high 46.44376 56.04234 50.91299 High range\n(340 - 390 Hz)
## 3  low 41.17250 48.98075 44.85125  Low range\n(220 - 330 Hz)
ga0 <- ggplot(pframe, aes(x=trt3, y=predMean))+
     geom_point()+ 
     labs( y = expression ("Sonication amplitude "(m~s^{-2})), x = "Frequency range for reward") + 
     geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1)+
     theme(axis.text.x = element_text(angle = 45, hjust = 1), 
           legend.position = 'none') +
     theme_classic() + 
  annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 58, label=c("a", "b", "c"),
                color="black") 
ga0

ggsave(plot = ga0, filename = file.path(figDir, "SonicationAmpPredsAndCI_unadjusted.pdf"), width = 5, height = 4)


# make plot using bonferroni adjustments

ga0 <- ggplot(pframe, aes(x=trt3, y=predMean))+
     geom_point()+ 
     labs( y = expression ("Sonication amplitude "(m~s^{-2})), x = "Frequency range for reward") + 
     geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1)+
     theme(axis.text.x = element_text(angle = 45, hjust = 1), 
           legend.position = 'none') +
     theme_classic() + 
  annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 58, label=c("a", "b", "a"),
                color="black") 
ga0

ggsave(plot = ga0, filename = file.path(figDir, "SonicationAmpPredsAndCI_BonfAdjusted.pdf"), width = 5, height = 4)

Visualize aggregated data for amplitude

aggdata <- aggregate(sl$amp_acc, by=list(colNum = sl$colNum,trialNum =sl$trialNum, trt = sl$trt), FUN=mean, na.rm=TRUE)
colnames(aggdata)[colnames(aggdata) == "x"] = "amp"

aggdata_sd <- aggregate(sl$amp_acc, by=list(colNum = sl$colNum,trialNum =sl$trialNum, trt = sl$trt), FUN=sd, na.rm=TRUE)
colnames(aggdata_sd)[colnames(aggdata_sd) == "x"] = "amp_sd"

aggdata = merge(aggdata, aggdata_sd)
aggdata = aggdata[order(aggdata$colNum, aggdata$trialNum, decreasing = FALSE), ]
agg_sm = aggdata[aggdata$trialNum <= 2, ]
rownames(agg_sm) = 1:nrow(agg_sm)
agg_sm
##                colNum trialNum    trt      amp    amp_sd
## 1              blue_4        1   full 40.03641 27.005004
## 2              gold_3        1   full 32.58012 13.395404
## 3           goldred_4        1   full 27.49328  8.575236
## 4             green_4        1   full 44.14520 45.654277
## 5             green_4        2 full_2 19.19985 21.543257
## 6              lime_5        1   full 72.53866 33.940021
## 7          limeblue_5        1   full 49.20279 24.837998
## 8          limeblue_5        2 full_2 59.04592 36.444346
## 9          limegold_5        1   full 34.95839 26.390545
## 10        limegreen_5        1   full 52.56061 32.949195
## 11       limeorange_5        1   full 43.28056 22.966005
## 12       limeorange_5        2 full_2 50.75113 25.907474
## 13         limepink_5        1   full 49.57493 19.711454
## 14       limepurple_5        1   full 43.98249 22.702402
## 15       limepurple_5        2    low 39.95722 19.359774
## 16 limepurpleyellow_5        1   full 34.54496 17.151702
## 17          limered_5        1   full 42.27701 20.963674
## 18          limered_5        2    low 40.63711 20.895896
## 19       limesilver_5        1   full 26.24910  7.336471
## 20        limewhite_5        1   full 37.65233 22.862123
## 21        limewhite_5        2 full_2 40.91118 19.936507
## 22       limeyellow_5        1   full 43.51579 23.156058
## 23           orange_3        1   full 32.25152 18.932874
## 24           orange_5        1   full 55.08235 20.197053
## 25       orangeblue_5        1   full 51.97426 23.301810
## 26      orangegreen_5        1   full 58.59866 33.410963
## 27       orangepink_5        1   full 50.98635 28.370137
## 28     orangepurple_5        1   full 34.53312 16.739467
## 29           purple_3        1   full 53.61466 21.935890
## 30          redblue_4        1   full 28.88506 11.312274
## 31          redblue_4        2 full_2 48.61796 25.277613
## 32         redgreen_5        1   full 78.72927 46.273733
## 33         redgreen_5        2 full_2 70.38830 47.080737
## 34          redpink_5        1   full 74.20294 37.536798
## 35          redpink_5        2    low 56.85081 31.951114
## 36        redpurple_5        1   full 58.06343 40.818616
## 37        redpurple_5        2   high 53.13637 26.813039
## 38           silver_5        1   full 46.34517 30.280907
## 39            white_4        1   full 47.76100 21.162920
## 40            white_4        2 full_2 46.45015 22.564532
## 41        whiteblue_5        1   full 66.11499 33.532056
## 42        whiteblue_5        2    low 67.00735 29.330027
## 43        whitegold_5        1   full 43.10847 26.148044
## 44        whitegold_5        2    low 45.20369 27.241996
## 45       whitegreen_4        1   full 38.87762 18.184245
## 46      whiteorange_5        1   full 61.63802 39.835205
## 47      whiteorange_5        2    low 45.85399 25.838373
## 48        whitepink_5        1   full 63.53760 35.030181
## 49        whitepink_5        2   high 80.19250 47.237615
## 50      whitepurple_5        1   full 47.77247 35.753931
## 51         whitered_5        1   full 48.37233 25.860700
## 52         whitered_5        2   high 74.79446 45.764818
## 53      whiteyellow_5        1   full 65.50425 32.808944
## 54      whiteyellow_5        2 full_2 56.00923 33.372943
## 55       yellowblue_5        1   full 52.45834 30.688932
## 56       yellowblue_5        2   high 79.36156 41.580725
## 57      yellowgreen_5        1   full 61.58710 33.835105
## 58      yellowgreen_5        2 full_2 79.18968 27.697874
## 59     yelloworange_5        1   full 73.32616 36.109449
## 60     yelloworange_5        2   high 54.48161 24.329362
## 61       yellowpink_5        1   full 52.74199 24.869065
## 62       yellowpink_5        2    low 57.76271 29.939244
## 63     yellowpurple_5        1   full 62.77375 29.437162
## 64     yellowpurple_5        2    low 68.98565 35.865129
## 65        yellowred_5        1   full 62.77866 39.753892
## 66        yellowred_5        2 full_2 47.97187 28.376661
agg_sm[duplicated(data.frame(agg_sm$colNum, agg_sm$trialNum)), ]
## [1] colNum   trialNum trt      amp      amp_sd  
## <0 rows> (or 0-length row.names)
ggplot(agg_sm, aes(x = trt, y = amp, fill = trialNum > 1)) +
  geom_boxplot(alpha = 0.2) + 
  geom_point(aes(color = trialNum>1)) + 
  geom_line(aes(group = colNum))

diffdf <- sapply(unique(agg_sm$colNum), FUN = function(x){
  tmp = agg_sm[agg_sm$colNum == x, ]
  if(nrow(tmp) <= 1)
    diff = NA
  else
    diff = tmp$amp[tmp$trialNum == 2] - tmp$amp[tmp$trialNum == 1]
  return(diff)
})

trtDF = sapply(unique(agg_sm$colNum), FUN = function(x){
  tmp = agg_sm[agg_sm$colNum == x, ]
  ttrs = paste(tmp$trt[tmp$trialNum == 1], tmp$trt[tmp$trialNum == 2], sep = "_")
  return(ttrs)
})

buzzdiffs = data.frame(trtDF, diffdf)

ggplot(buzzdiffs, aes(x = trtDF, y= diffdf)) + 
  geom_boxplot() + 
  geom_point()+ 
  labs(y = "amplitude difference m/s/s")
## Warning: Removed 20 rows containing non-finite values (stat_boxplot).
## Warning: Removed 20 rows containing missing values (geom_point).

agg2 = aggregate(sl$amp_acc, by=list(colNum = sl$colNum, fullTrt = sl$trt == "full", trt = sl$trt), FUN=mean, na.rm=TRUE)
colnames(agg2)[colnames(agg2) == "x"] = "amp"
agg2$trt = as.character(agg2$trt)

diffdf <- t(as.data.frame(t(sapply(unique(agg2$colNum), FUN = function(x){
  tmp = agg2[agg2$colNum == x, ]
  if(nrow(tmp) <= 1)
    return(NA)
  if (length(unique(tmp$trt)) > 2){
    tmp = tmp[tmp$trt != "full_2", ]
  }
  diff =  tmp$amp[!tmp$fullTrt] - tmp$amp[tmp$fullTrt]
  return(diff)

}))))

length(diffdf)
## [1] 43
trtDF = sapply(unique(agg2$colNum), FUN = function(x){
  tmp = agg2[agg2$colNum == x, ]
  if (length(unique(tmp$trt)) > 2){
    tmp = tmp[tmp$trt != "full_2", ]
  }
  ttrs = paste(tmp$trt[tmp$fullTrt], tmp$trt[!tmp$fullTrt], sep = "_")
  return(ttrs)
})

length(trtDF)
## [1] 43
buzzdiffs = data.frame(trtDF, diffdf)
tapply(buzzdiffs$diffdf, INDEX = buzzdiffs$trtDF, mean)
##       full_ full_full_2   full_high    full_low 
##          NA    2.871714    9.194234    3.249068
ggplot(droplevels(buzzdiffs[buzzdiffs$trtDF != "full_", ]), aes(x = trtDF, y= as.numeric(diffdf))) + 
  geom_boxplot() + 
  geom_point()+ 
  labs(y = "amplitude difference m/s/s")

Plot amplitude vs. frequency for initial trial, accounting for size

ggplot(sl[sl$trialNum == 1, ], aes(x = freq, y = amp_acc)) + 
  geom_point(position = position_jitter(width = 3), alpha = 0.2) + 
  theme_classic() + 
  geom_smooth(aes(group = trt2)) + 
  labs( y = expression ("Sonication amplitude "(m~s^{-2})), x = "Frequency (Hz)")
## `geom_smooth()` using method = 'gam'

centerFreq = scale(sl[sl$trialNum == 1, "freq"])

m1 <- lmer(log(amp_acc)~centerFreq + I(centerFreq^2) + I(centerFreq^3) + IT_imputed + (1|beeCol), data = sl[sl$trialNum == 1, ])
summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(amp_acc) ~ centerFreq + I(centerFreq^2) + I(centerFreq^3) +  
##     IT_imputed + (1 | beeCol)
##    Data: sl[sl$trialNum == 1, ]
## 
## REML criterion at convergence: 3526.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5361 -0.6344  0.0693  0.6857  2.4239 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeCol   (Intercept) 0.08395  0.2897  
##  Residual             0.32877  0.5734  
## Number of obs: 1971, groups:  beeCol, 42
## 
## Fixed effects:
##                 Estimate Std. Error t value
## (Intercept)      2.00589    0.57016   3.518
## centerFreq       0.37124    0.03004  12.357
## I(centerFreq^2) -0.08590    0.01227  -7.003
## I(centerFreq^3) -0.03835    0.01055  -3.636
## IT_imputed       0.43530    0.13776   3.160
## 
## Correlation of Fixed Effects:
##             (Intr) cntrFr I(F^2) I(F^3)
## centerFreq  -0.047                     
## I(cntrFr^2) -0.093  0.081              
## I(cntrFr^3) -0.018 -0.844  0.065       
## IT_imputed  -0.996  0.043  0.071  0.019
plot(m1)

# predict -- using mean IT span
ppdf <- data.frame(centerFreq = sort(unique(centerFreq)), IT_imputed = ITmean, beeCol = 99999, acc_amp = 0)
# exponentiate to get back to original scale
ppdf$predAmp = exp(predict(m1, newdata = ppdf, type = "response", re.form = NA))



cf_Unscaled = ppdf$centerFreq * attr(centerFreq, 'scaled:scale') + attr(centerFreq, 'scaled:center')
ggplot(ppdf, aes(x = cf_Unscaled, y = predAmp))+ 
  geom_line() + 
  labs(x = "Sonication Frequency (Hz)", y = expression ("Sonication amplitude "(m~s^{-2})) ) + 
  geom_point(data = sl[sl$trialNum == 1, ], 
             aes(x = freq, y = amp_acc), 
             alpha = 0.2, position = position_jitter(width =2), pch =16, size = .5) 

ggsave(filename = file.path(figDir, "FreqVsAmp_1stTrial_rawData.pdf"), width = 4, height = 3)

g22 <- ggplot(ppdf, aes(x = cf_Unscaled, y = predAmp))+ 
  geom_line() + 
  labs(x = "Sonication Frequency (Hz)", y = expression ("Predicted sonication amplitude "(m~s^{-2})) ) 
g22

ggsave(filename = file.path(figDir, "FreqVsAmp_1stTrial.pdf"), width = 4, height = 3)
sessionInfo()
## R version 3.4.2 (2017-09-28)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.2
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] viridis_0.4.0     viridisLite_0.2.0 effects_4.0-0    
##  [4] carData_3.0-0     plyr_1.8.4        multcomp_1.4-8   
##  [7] TH.data_1.0-8     MASS_7.3-47       survival_2.41-3  
## [10] mvtnorm_1.0-6     sjPlot_2.4.0      lme4_1.1-14      
## [13] Matrix_1.2-11     reshape2_1.4.2    ggplot2_2.2.1    
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-131       pbkrtest_0.4-7     RColorBrewer_1.1-2
##  [4] rprojroot_1.2      tools_3.4.2        TMB_1.7.11        
##  [7] backports_1.1.1    R6_2.2.2           sjlabelled_1.0.5  
## [10] DT_0.2             mgcv_1.8-20        lazyeval_0.2.1    
## [13] colorspace_1.3-2   nnet_7.3-12        tidyselect_0.2.3  
## [16] gridExtra_2.3      mnormt_1.5-5       compiler_3.4.2    
## [19] sandwich_2.4-0     labeling_0.3       scales_0.5.0      
## [22] lmtest_0.9-35      psych_1.7.8        blme_1.0-4        
## [25] stringr_1.2.0      digest_0.6.12      foreign_0.8-69    
## [28] minqa_1.2.4        rmarkdown_1.8      stringdist_0.9.4.6
## [31] pkgconfig_2.0.1    htmltools_0.3.6    pwr_1.2-1         
## [34] htmlwidgets_0.9    rlang_0.1.4        shiny_1.0.5       
## [37] bindr_0.1          zoo_1.8-0          dplyr_0.7.4       
## [40] magrittr_1.5       modeltools_0.2-21  bayesplot_1.4.0   
## [43] Rcpp_0.12.13       munsell_0.4.3      abind_1.4-5       
## [46] prediction_0.2.0   stringi_1.1.6      yaml_2.1.14       
## [49] merTools_0.3.0     snakecase_0.5.1    grid_3.4.2        
## [52] parallel_3.4.2     sjmisc_2.6.2       forcats_0.2.0     
## [55] lattice_0.20-35    ggeffects_0.2.2    haven_1.1.0       
## [58] splines_3.4.2      sjstats_0.12.0     knitr_1.17        
## [61] codetools_0.2-15   stats4_3.4.2       glue_1.2.0        
## [64] evaluate_0.10.1    modelr_0.1.1       nloptr_1.0.4      
## [67] httpuv_1.3.5       gtable_0.2.0       purrr_0.2.4       
## [70] tidyr_0.7.2        assertthat_0.2.0   mime_0.5          
## [73] coin_1.2-1         xtable_1.8-2       broom_0.4.2       
## [76] survey_3.32-1      coda_0.19-1        tibble_1.3.4      
## [79] arm_1.9-3          glmmTMB_0.1.4      bindrcpp_0.2